1 Area of Interest

2 Landsat Tiles

The dimensions of the stacked image are 7,811 rows, 7,861 columns, 59,996,291 cells and 6 layers. The CRS is UTM on the WGS84 datum. The resolution is 30m x 30m.

The dimensions of the stacked image are 340 rows, 346 columns, 117,640 cells and 6 layers. The CRS is UTM on the WGS84 datum. The resolution is 30m x 30m.

3 Plots

3.1 Non-Contrasted

3.1.1 Natural Color

3.1.2 Color Infared

3.1.3 False Color Water Focus

3.1.4 Land From Water Focus

3.2 Contrasted

3.2.1 Natural Color

3.2.2 Color Infared

3.2.3 False Color Water Focus

3.2.4 Land From Water Focus

4 Flooding

4.1 Raster Algebra

A we can see, the NWDI, MNDWI, and WRI are all very similar in the way they pick up water cells with the same color, but they deviate differently in the dry ground cells. NDVI picks up water in a different color, and land in the opposite as well. The SWI seems to pick up only water cells and no land cells.

4.2 Thresholding

SWI has multiple NA values so we will have to replace those.

5 Clustering

This shoes us they are extracted as a matrix of size 117640 x 6

##    
##         1     2     3     4     5     6
##   0   488 17474  8361 53172 22289   655
##   1 13387     0    50   108  1646    10

6 Flooding Stats

Flood Area for Different Methods
Flood Area m^2
NDVI 5999400
NDWI 6490800
MNDWI 10745100
WRI 7622100
SWI 13680900
K.Means 12487500

6.1 Flood Uncertainty

6.1.1 Base R

6.1.2 Open Steet Map

Some of the values in the raster are not even because there may be 1, 3, or 5 raster layers which picked up areas of the map as being flooded, which would result in the sum of layers not being even.

7 Appendix

knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE)
library(tidyverse)
library(sf)
library(raster)
library(getlandsat)
library(mapview)
library(osmdata)
library(knitr)
bb = read_csv('../data/uscities.csv') %>%
  filter(city == 'Palo') %>%
  st_as_sf(coords = c('lng','lat'), crs = 4326) %>%
  st_transform(5070) %>%
  st_buffer(5000) %>%
  st_bbox() %>%
  st_as_sfc()
meta = read_csv('../data/palo-flood.csv')

files = lsat_scene_files(meta$download_url) %>%
  filter(grepl(paste0('B',1:6,'.TIF$', collapse = '|'), file)) %>%
  arrange(file) %>% 
  pull(file)

st = sapply(files, lsat_image) 

b = stack(st) %>% setNames(paste0('band', 1:6))
cropper = bb %>% st_as_sf() %>% st_transform(crs(b))

cr_b = crop(b, cropper)

r <- cr_b %>% setNames(c('coastal', 'blue', 'green', 'red', 'nir', 'swir'))

plotRGB(r, r = 4, g= 3, b = 2)
plotRGB(r, r = 5, g= 4, b = 2)
plotRGB(r, r = 5, g= 6, b = 4)
plotRGB(r, r = 7, g= 5, b = 3)
plotRGB(r, r = 4, g= 3, b = 2, stretch = 'hist')
plotRGB(r, r = 5, g= 4, b = 2, stretch = 'hist')
plotRGB(r, r = 5, g= 6, b = 4, stretch = 'lin')
plotRGB(r, r = 7, g= 5, b = 3, stretch = 'lin')
ndvi = (r$nir - r$red)/(r$nir + r$red)

ndwi = (r$green - r$nir)/(r$green + r$nir)

mndwi = (r$green - r$swir)/(r$green + r$swir)

wri = (r$green + r$red)/(r$nir + r$swir)

swi = 1/(sqrt(r$blue-r$swir))

rs <- stack(ndvi, ndwi, mndwi, wri, swi) %>% setNames(c('NDVI', 'NDWI', 'MNDWI', 'WRI', 'SWI'))

plot(rs, col = colorRampPalette(c("blue", "white", "red"))(256))

ndvi_t = function(x) {ifelse(x<=0,1,0)}
ndwi_t = function(x) {ifelse(x>=0,1,0)}
mndwi_t = function(x) {ifelse(x>=0,1,0)}
wri_t = function(x) {ifelse(x>=1,1,0)}
swi_t = function(x) {ifelse(x<=5,1,0)}


ndvi_f = calc(ndvi, ndvi_t)
ndwi_f = calc(ndwi, ndwi_t)
mndwi_f = calc(mndwi, mndwi_t)
wri_f = calc(wri, wri_t)
swi_f = calc(swi, swi_t)

f_stack = stack(c(ndvi_f, ndwi_f, mndwi_f, wri_f, swi_f)) %>% setNames(c('NDVI', 'NDWI', 'MNDWI', 'WRI', 'SWI'))
plot(f_stack, col = colorRampPalette(c("white","blue"))(256))
#sum(is.na(values(ndvi_f)))
#sum(is.na(values(ndwi_f)))
#sum(is.na(values(mndwi_f)))
#sum(is.na(values(wri_f)))
#sum(is.na(values(swi_f)))
values(swi_f) <- ifelse(is.na(values(swi_f)), 0, values(swi_f) )

#sum(is.na(values(swi_f)))
vals <- getValues(r)

#dim(vals)
vals <- na.omit(vals)
vs <- scale(vals)
idx <- 1:dim(vals)[1]
set.seed(20200905)

rast.clust <- kmeans(vs, 6, iter.max = 100)

kmeans_raster <- rs$MNDWI
values(kmeans_raster) <- rast.clust$cluster

plot(kmeans_raster)
table(getValues(swi_f), getValues(kmeans_raster))
k_thresh = function(x) {ifelse(x==1,1,0)}

k_flood = calc(kmeans_raster, k_thresh) %>% setNames('K-Means')

f_stack <- addLayer(f_stack, k_flood)
plot(f_stack, col = colorRampPalette(c("white","blue"))(256))
flood_area <- cellStats(f_stack, 'sum') * (30^2)
flood_area_km <- flood_area*0.000001

#names(flood_area)

kable(data.frame(flood_area), col.names = 'Flood Area m^2', caption = 'Flood Area for Different Methods')
rast_uncert <- calc(f_stack, function(x){sum(x)}) %>% setNames('Flood Uncertainty')
plot(rast_uncert, col = blues9, main = 'Flood Uncertainty')
values(rast_uncert) <- ifelse(values(rast_uncert)==0, NA, values(rast_uncert))
mapview(rast_uncert)